clear; clc; close all;

% 参数设置
a = 0.1;
b = 0.1;
c = 0.2;         % Fig.8(b) 中 c=0.2
k = 1.85;
x0_fixed = 0.1;  % 固定 x0

q0_min = -20;
q0_max = 20;
num_q0 = 500;
q0_list = linspace(q0_min, q0_max, num_q0);

% 迭代参数
transient = 1000;  % 舍去初始暂态步数
N_plot = 900;      % 用于绘制分岔图的记录步数
N_LE = 1000;       % 用于计算 Lyapunov 指数的有效迭代步数

% 用于存储每个 q0 下记录的 x 与 q 值
bifurcationX = cell(1, num_q0);
bifurcationQ = cell(1, num_q0);

% 用于存储每个 q0 对应的 LEs 指数
LE1_vals = zeros(1, num_q0);
LE2_vals = zeros(1, num_q0);

% 对每个初始 q0 进行迭代
for idx = 1:num_q0
    q0 = q0_list(idx);
    % 初始化状态
    x = x0_fixed;
    q = q0;
    
    % 先进行 transient 步的迭代，去除暂态影响
    for t = 1:transient
        [x, q] = mptm(x, q, a, b, c, k);
    end
    
    % 记录后续 N_plot 步的 x 与 q 值
    x_values = zeros(1, N_plot);
    q_values = zeros(1, N_plot);
    for t = 1:N_plot
        [x, q] = mptm(x, q, a, b, c, k);
        x_values(t) = x;
        q_values(t) = q;
    end
    bifurcationX{idx} = x_values;
    bifurcationQ{idx} = q_values;
    
    % 计算该初始条件下的 LEs 指数
    [LE1, LE2] = LEs(@mptm, x0_fixed, q0, a, b, c, k, N_LE, transient);
    LE1_vals(idx) = LE1;
    LE2_vals(idx) = LE2;
end

% 上图：LEs 指数随 q0 的变化
subplot(2,1,1);
hold on; grid on;
plot(q0_list, LE1_vals, 'b', 'LineWidth', 1.5);
plot(q0_list, LE2_vals, 'r', 'LineWidth', 1.5);
xlabel('q_0');
ylabel('LEs 指数');
legend('LE1','LE2','Location','Best');
title('LEs 指数 (c=0.2, k=1.85)');

% 下图：分岔图
subplot(2,1,2);
hold on; grid on;
for idx = 1:num_q0
    q_val = q0_list(idx);
    x_data = bifurcationX{idx};
    q_data = bifurcationQ{idx};
    
    % 绘制 x 值（红色）
    plot(q_val * ones(size(x_data)), x_data, '.r', 'MarkerSize', 1);
    % 绘制 q 值（蓝色）
    plot(q_val * ones(size(q_data)), q_data, '.b', 'MarkerSize', 1);
end
xlabel('q_0');
ylabel('x (红色) 与 q (蓝色)');
title('分岔图 (c=0.2, k=1.85)');
